{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第一周 线性方程组与矩阵\n",
    "- 陈怀琛电子邮址是：hchchen1934@vip.163.com\n",
    "- 杨威电子邮址是：weiyang@mail.xidian.edu.cn\n",
    "\n",
    "### 线性方程应用\n",
    "应用：某食品是由原材料甲、乙、丙、丁四种食品配置成的,\n",
    "问四种原材料占比多少能配出某食品？\n",
    "\n",
    "|成分|甲|乙|丙|丁|某食品|\n",
    "|-|-|-|-|-|-|\n",
    "|蛋白质(%)|20|16|10|15|15|\n",
    "|脂肪(%)|3|8|2|5|5|\n",
    "|碳水化合物(%)|10|25|20|5|12|\n",
    "|未知数|$x_1$|$x_2$|$x_3$|$x_4$|\n",
    "\n",
    "$\\begin{cases}\n",
    "x_1 + x_2 + x_3 + x_4 = 1\\\\\n",
    "20x_1 + 16x_2 + 10x_3 + 15x_4 = 15\\\\\n",
    "3x_1 + 8x_2 + 2x_3 + 5x_4 = 5\\\\\n",
    "10x_1 + 25x_2 + 20x_3 + 5x_4 = 12\n",
    "\\end{cases}$\n",
    "\n",
    "$A=\\begin{bmatrix}1,1,1,1\\\\20,16,10,15\\\\3,8,2,5\\\\10,25,20,5\\end{bmatrix}$ A是系数矩阵\n",
    "\n",
    "$B=\\begin{bmatrix}1\\\\15\\\\5\\\\12\\end{bmatrix}$ B是常数项\n",
    "\n",
    "$C[A,B]=\\begin{bmatrix}1,1,1,1,1\\\\20,16,10,15,15\\\\3,8,2,5,5\\\\10,25,20,5,12\\end{bmatrix}$ C是增广矩阵\n",
    "\n",
    "> 1. 第一个方程四种食品配比占比100%等于1\n",
    "> 2. 二、三、四左右都乘以100\n",
    "> 3. 这些都是约束条件 \n",
    "### 线性方程应用Numpy求解：\n",
    "AX = B 可以写成\n",
    "\n",
    "$\\begin{bmatrix}1,1,1,1\\\\20,16,10,15\\\\3,8,2,5\\\\10,25,20,5\\end{bmatrix}$ ×\n",
    "$\\begin{bmatrix}x_1\\\\x_2\\\\x_3\\\\x_4\\end{bmatrix}$ = \n",
    "$\\begin{bmatrix}1\\\\15\\\\5\\\\12\\end{bmatrix}$\n",
    "\n",
    "**下面是用numpy求解x的过程** \\\n",
    "参考资料：https://blog.csdn.net/wqqgo/article/details/85101888"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "变量x的集合：[0.10306748 0.21472393 0.14601227 0.53619632]\n",
      "常数项 b 检测：[ 1. 15.  5. 12.]\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "a = np.array([[1,1,1,1],[20,16,10,15],[3,8,2,5],[10,25,20,5]])  # 系数矩阵\n",
    "b = np.array([1,15,5,12])  # 常数项 \n",
    "# a+b是 增广矩阵\n",
    "x = np.linalg.solve(a, b)  # 求解线性方程组 高斯消元后x1，x2，x3分别是 Ax = b => x = b / A\n",
    "check = np.dot(a,x)  # 矩阵的点积Numpy方法 check a * x 是不是 等于 b\n",
    "print(f\"变量x的集合：{x}\")\n",
    "print(f\"常数项 b 检测：{check}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 线性代数的任务之一\n",
    "#### 线性方程组解的情况：\n",
    "- 有解\n",
    "    - 有多解 -> 解集的性质\n",
    "    - 有唯一解 -> 合理解、不合理解\n",
    "- 无解 -> 找出近似值\n",
    "\n",
    "#### 线性方程组的解有三中类型：\n",
    "1. 适定方程组：存在着唯一的一组解\n",
    "2. 欠定方程组：其解存在但不唯一\n",
    "3. 超定方程组：不存在精确解，可以求出其近似值\n",
    "\n",
    "#### 方程组的几何意义是焦点也就是斜率\n",
    "1. 两个未知数的方程就是集合意义上的焦点\n",
    "2. 两个未知数有三个方程，发现两两相交只能求出超定方程近似解\n",
    "3. 三个未知数意义就是三维平面的相交，有可能有一个点，有可能是无数个三个平面交与一条直线 \n",
    "\n",
    "#### 高斯消元法与阶梯型方程组\n",
    "1. 互换两个方程的位置\n",
    "2. 用一个非零数k乘某个方程；称为数乘变换\n",
    "3. 把一个方程的k倍加到另外一个方程上。称为消元变换\n",
    "    - 线性方程组的同解变换\n",
    "\n",
    "#### 利用矩阵最简形解方程组\n",
    "##### 矩阵-行最简形：\n",
    "1. 阶梯矩阵 row 每个阶梯1个阶梯,\n",
    "\n",
    "$\\begin{bmatrix}1,2,3,4,5\\\\0,2,3,7,9\\\\0,0,0,3,8\\\\0,0,0,0,0\\end{bmatrix}$ 是阶梯矩阵不是最简形\n",
    "2. 行最简形需要满足3个条件\n",
    "    - 是阶梯矩阵\n",
    "    - 每一个row第一个数字是1，前面没有或者是0\n",
    "    - 这个1所在列除了它都是0（列首元）\n",
    "\n",
    "$\\begin{bmatrix}1,0,0,4,5,7\\\\0,1,0,7,9,3\\\\0,0,1,0,3,8\\end{bmatrix}$\n",
    "3. 矩阵的最简形和原增广矩阵是同解\n",
    "\n",
    "[矩阵最简单形状-wiki](https://zh.wikipedia.org/zh-hans/%E9%98%B6%E6%A2%AF%E5%BD%A2%E7%9F%A9%E9%98%B5)\n",
    "\n",
    "##### 矩阵最简形具体功能：\n",
    "1. 解线性方程（系数矩阵、常数项）\n",
    "2. 求矩阵的秩（zhi）；秩含义：化为最简形后有几row非全部零，就有几个约束条件（简单解释）\n",
    "3. 求矩阵行最简形首元所在的列数，第一列1,就是$x_1$第三列1就是$x_3$\n",
    "\n",
    "$\\begin{bmatrix}1,0,0,4,5，7\\\\0,1,0,7,9,3\\\\0,0,1,0,3,8\\\\0,0,0,0,0,0\\end{bmatrix}$ 秩是3个\n",
    "\n",
    "##### 应用计算方程组案列1：\n",
    "$\\begin{cases}\n",
    "2x_1 - 2x_2 + 2x_3 + 6x_4 = -16\\\\\n",
    "2x_1 - x_2 + 2x_3 + 4x_4 = -10\\\\\n",
    "3x_1 - x_2 + 4x_3 + 4x_4 = -11\\\\\n",
    "x_1 + x_2 - x_3 + 3x_4 = -12\n",
    "\\end{cases}$\n",
    "结果：\n",
    "$\\begin{bmatrix}\n",
    "1,0,0,0,11\\\\\n",
    "0,1,0,0,-8\\\\\n",
    "0,0,1,0,-6\\\\\n",
    "0,0,0,1,-7\n",
    "\\end{bmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x变量： [11. -8. -6. -7.]\n",
      "[-16. -10. -11. -12.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "# a * x = b  矩阵的点积？为什么这样乘可以和上文x1、x2、x3、x4 对应？\n",
    "a = np.array([[2,-2,2,6],[2,-1,2,4],[3,-1,4,4],[1,1,-1,3]])  # 系数矩阵\n",
    "b = np.array([-16,-10,-11,-12])  # 常数项\n",
    "x = np.linalg.solve(a,b)\n",
    "print(\"x变量：\",x)\n",
    "check = np.dot(a,x)  # 矩阵的点积Numpy方法\n",
    "print(check)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 应用计算方程组案列2：\n",
    "设方程组的系数矩阵A，b如下\n",
    "$A=\\begin{bmatrix}\n",
    "-2&-2&2&2&-2\\\\\n",
    " 1&-5&1&-3&-1\\\\\n",
    "-1&2&-5&6&5\\\\\n",
    "-1&2&1& 0 &-1\n",
    "\\end{bmatrix}$\n",
    "，$b=\\begin{bmatrix}-2,-1,2,0\\end{bmatrix}$\n",
    "\n",
    "得到解：\n",
    "$c=\\begin{bmatrix}\n",
    "1&0&0&0&0&-2/9\\\\\n",
    "0&1&0&0&0&2/9\\\\\n",
    "0&0&1&-1&0&-2/3\\\\\n",
    "0&0&0&0&1&-1/3\n",
    "\\end{bmatrix}$\n",
    "\n",
    "1. ip = 1 2 3 5  # 只有第一列、二列、三列、五列是行首元 \n",
    "2. 秩有四个行非零，五个未知数（列）\n",
    "3. 第三行有一个约束条件两个未知数 $x_3,x_4$\n",
    "\n",
    "$\\begin{cases}\n",
    "x_1 = -2/9\\\\\n",
    "x_2 = 2/9\\\\\n",
    "x_5 = -1/3\\\\\n",
    "\\end{cases}$\n",
    "且 \n",
    "$\\begin{cases}\n",
    "x_3 + x_4 = -2/3\n",
    "\\end{cases}$\n",
    "第三行无穷解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "LinAlgError",
     "evalue": "Singular matrix",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mLinAlgError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-3-06c8b1e9c68b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m160\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m40\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m210\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m330\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mb\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      5\u001b[0m \u001b[1;31m# Singular matrix 是一个欠定方程 解有无穷多个，用下面自己编写的函数求解\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\numpy\\linalg\\linalg.py\u001b[0m in \u001b[0;36msolve\u001b[1;34m(a, b)\u001b[0m\n\u001b[0;32m    401\u001b[0m     \u001b[0msignature\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'DD->D'\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32melse\u001b[0m \u001b[1;34m'dd->d'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    402\u001b[0m     \u001b[0mextobj\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 403\u001b[1;33m     \u001b[0mr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgufunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    404\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    405\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\numpy\\linalg\\linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[1;34m(err, flag)\u001b[0m\n\u001b[0;32m     95\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     96\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 97\u001b[1;33m     \u001b[1;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Singular matrix\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     98\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     99\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mLinAlgError\u001b[0m: Singular matrix"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "A = np.array([[1,-1,0,0],[0,1,-1,0],[0,0,1,-1],[-1,0,0,1]])\n",
    "b = np.array([160,-40,210,-330])\n",
    "x = np.linalg.solve(A,b)\n",
    "# Singular matrix 是一个欠定方程 解有无穷多个，用下面自己编写的函数求解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.   0.   0.  -1. 330.]\n",
      " [  0.   1.   0.  -1. 170.]\n",
      " [  0.   0.   1.  -1. 210.]\n",
      " [  0.   0.   0.   0.   0.]]\n"
     ]
    }
   ],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Created on Fri Mar 15 15:05:53 2019\n",
    "求解行最简行公式\n",
    "@author: MrZjjPolarBear\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "def rref(arbmat):  # 效仿MatLab函数 Reduced Row Echelon Form\n",
    "    \"\"\" Convert an arbitrary matrix to a simplest matrix \"\"\"\n",
    "    arbmat = arbmat.astype(float)\n",
    "    row_number, column_number = arbmat.shape\n",
    "    if row_number == 1:\n",
    "        if arbmat[0, 0] != 0:\n",
    "            return (arbmat/arbmat[0, 0])\n",
    "        else:\n",
    "            return arbmat\n",
    "    else:\n",
    "        rc_number = min(row_number, column_number)\n",
    "        anarbmat = arbmat.copy()\n",
    "        r = 0\n",
    "        for n in range(rc_number):\n",
    "            s_row = -1\n",
    "            for i in arbmat[r:row_number, n]:\n",
    "                s_row += 1\n",
    "                if abs(i) > 1e-10:\n",
    "                    anarbmat[r, :] = arbmat[s_row+r, :]\n",
    "                    for j in range(r, row_number):\n",
    "                        if j < s_row+r:\n",
    "                            anarbmat[j+1, :] = arbmat[j, :]\n",
    "                    arbmat = anarbmat.copy()\n",
    "            if abs(anarbmat[r, n]) > 1e-10:\n",
    "                anarbmat[r, :] = anarbmat[r, :] / anarbmat[r, n]\n",
    "                for i in range(row_number):\n",
    "                    if i != r:\n",
    "                        anarbmat[i, :] -= anarbmat[i, n]*anarbmat[r, :]\n",
    "            arbmat = anarbmat.copy()\n",
    "            if abs(arbmat[r, n]) < 1e-10:\n",
    "                r = r\n",
    "            else:\n",
    "                r = r + 1\n",
    "        for m in range(column_number):\n",
    "            if abs(arbmat[-1, m]) > 1e-10:\n",
    "                arbmat[-1, :] = arbmat[-1, :]/arbmat[-1, m]\n",
    "                for i in range(row_number-1):\n",
    "                    arbmat[i, :] -= arbmat[i, m]*arbmat[-1, :]\n",
    "                break\n",
    "            \n",
    "        return arbmat\n",
    "a = np.matrix([[1,-1,0,0,160],[0,1,-1,0,-40],[0,0,1,-1,210],[-1,0,0,1,-330]])\n",
    "# 160,-40,210,-330 是常数项和MATLAB rref(A,b) -> rsmat([A,b])\n",
    "a1 = rref(a)\n",
    "print(a1)  # 行最简形"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 随机生产m=4，n=3，取值范围0-10 矩阵 且求矩阵的秩"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "随机生产的矩阵是： \n",
      " [[0 0 0]\n",
      " [0 0 1]\n",
      " [0 1 1]\n",
      " [0 1 1]]\n",
      "此矩阵的约束条件秩是: 2\n",
      "行最简性是： \n",
      " [[0. 1. 0.]\n",
      " [0. 0. 1.]\n",
      " [0. 0. 0.]\n",
      " [0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "b = np.random.randint(0,2,(4,3))\n",
    "c = np.linalg.matrix_rank(b)\n",
    "print(\"随机生产的矩阵是：\",\"\\n\",b)\n",
    "print(\"此矩阵的约束条件秩是:\",c)  # rank秩\n",
    "a2 = rref(b)\n",
    "print(\"行最简性是：\",\"\\n\",a2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.        ,  0.        ,  1.33333333,  0.        ],\n",
       "       [ 0.        ,  1.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [-0.        , -0.        ,  1.        , -0.33333333, -0.        ]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([[1,1,1,1,0],[2,1,-1,3,0],[1,-2,1,1,0]])\n",
    "rref(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 练习题1.72 1.4题(d) 体现了的高斯消元算法思路经典啊！！！！\n",
    "\n",
    "$\n",
    "(d)=\\begin{bmatrix}\n",
    "1&1&1&1&0 \\\\\n",
    "2&1&-1&3&0  \\\\\n",
    "1&-2&1&1&0  \n",
    "\\end{bmatrix}\n",
    "$\n",
    "\n",
    "解题思路：\n",
    "1. 先从第一行为基准，从上而下、从左至右消元\n",
    "    - C1(2,:) = C1(2,:) - C1(2,1) / C1(1,1) * C1(1,:)  # 数字2代表第二行\n",
    "    - C1(3,:) = C1(3,:) - C1(3,1) / C1(1,1) * C1(1,:) \n",
    "    - C1(3,:) = C1(3,:) - C1(3,2) / C1(2,2) * C1(2,:) \n",
    "\n",
    "$C1=\\begin{bmatrix}\n",
    "1&1&1&1&0 \\\\\n",
    "2&1&-1&3&0  \\\\\n",
    "1&-2&1&1&0 \n",
    "\\end{bmatrix}$ 变成 \n",
    "$C=C1=\\begin{bmatrix}\n",
    "1&1&1&1&0 \\\\\n",
    "0&-1&-3&1&0  \\\\\n",
    "0&0&9&-3&0  \n",
    "\\end{bmatrix}$\n",
    "\n",
    "2. 发现是阶梯矩阵但不是行最简性，首元非1且列非除零外只有1\n",
    "3. 先从最后一行为基准，从下而上、从右至左消元\n",
    "    - C1(2,:) = C1(2,:) - C1(2,3) / C1(3,3) * C1(3,:) \n",
    "    - C1(1,:) = C1(1,:) - C1(1,3) / C1(3,3) * C1(3,:) \n",
    "    - C1(1,:) = C1(1,:) - C1(1,2) / C1(2,2) * C1(2,:) \n",
    "\n",
    "$C1=\\begin{bmatrix}\n",
    "1&0&0&1.33333&0 \\\\\n",
    "0&1&0&0&0  \\\\\n",
    "0&0&1&-0.3333&0  \n",
    "\\end{bmatrix}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 第一周总结：\n",
    "1. 了解二阶、三阶在笛卡尔坐标系图形及几何意义，欠定方程（无数解、约束条件不够）、适定方程（唯一解）、超定方程（无解、不存在精确解、近似解）\n",
    "2. 线性方程 - 》高斯消元 -》矩阵变换（行最简形） 初等行变换\n",
    "3. 行最简形：是阶梯矩阵且每一个row第一个是1且每一列只有1（阶梯上）\n",
    "4. 名次：矩阵的秩就是初等行变换后真正的约束条件\n",
    "5. 系数矩阵+常数项 = 增广矩阵\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第二周 矩阵的四则运算\n",
    "### 用初等行变换求逆矩阵\n",
    "[重要信息可以看到初等变换求逆矩阵全过程](https://www.hahack.com/math/math-matrix/)\n",
    "\n",
    "(AI)−→−−−初等变换(IA<sup>−1</sup>)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a是矩阵：\n",
      "[[ 1  3 -2]\n",
      " [-3 -6  5]\n",
      " [ 1  1 -1]]\n",
      "a的逆矩阵通过初等变换：\n",
      "[[1. 1. 3.]\n",
      " [2. 1. 1.]\n",
      " [3. 2. 3.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "a = np.matrix('1 3 -2; -3 -6 5; 1 1 -1')\n",
    "print(f\"a是矩阵：\\n{a}\")\n",
    "print(f\"a的逆矩阵通过初等变换：\\n{a.I}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 矩阵的加减、数乘、相乘（矢量积）、数量积、矩阵的逆矩matrix、array语法     \n",
    "#### matrix语法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "矩阵A 加 矩阵B: \n",
      " [[ 330  600  110]\n",
      " [ 550 1050  220]\n",
      " [ 650 1550  270]\n",
      " [ 950 1700  180]]\n",
      "矩阵A数乘 加 矩阵B数乘：\n",
      " [[177. 340.  59.]\n",
      " [295. 545. 118.]\n",
      " [345. 835. 147.]\n",
      " [495. 850.  98.]]\n",
      "矩阵A乘矩阵b：\n",
      "[[273]\n",
      " [339]\n",
      " [217]]\n",
      "矩阵A与矩阵c每个元素相乘：\n",
      "[[400  16   4  49]\n",
      " [100 324  25  36]\n",
      " [ 25  49 256   9]]\n",
      "A矩阵的逆矩阵是：\n",
      "[[ 3. -2.]\n",
      " [-1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "A = np.matrix([[150,250,50],[250,500,100],[300,700,120],[450,850,80]])\n",
    "B = np.matrix([[180,350,60],[300,550,120],[350,850,150],[500,850,100]])\n",
    "C = A + B\n",
    "print(f'矩阵A 加 矩阵B: \\n {C}')\n",
    "D = 0.1 * A + 0.9 * B\n",
    "print(f'矩阵A数乘 加 矩阵B数乘：\\n {D}')\n",
    "A = np.matrix([[20,4,2,7],[10,18,5,6],[5,7,16,3]])\n",
    "b = np.matrix([[8],[10],[5],[9]])  # 等价与 np.matrix('8;10;5;9')\n",
    "C = A * b\n",
    "print(f'矩阵A乘矩阵b：\\n{C}')\n",
    "c = np.matrix([[20,4,2,7],[10,18,5,6],[5,7,16,3]])\n",
    "D = np.multiply(A,c)\n",
    "print(f'矩阵A与矩阵c每个元素相乘：\\n{D}')\n",
    "A = np.matrix([[1,2],[1,3]])\n",
    "X = A.I\n",
    "print(f'A矩阵的逆矩阵是：\\n{X}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### array语法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "矩阵A 加 矩阵B: \n",
      " [[ 330  600  110]\n",
      " [ 550 1050  220]\n",
      " [ 650 1550  270]\n",
      " [ 950 1700  180]]\n",
      "矩阵A数乘 加 矩阵B数乘：\n",
      " [[177. 340.  59.]\n",
      " [295. 545. 118.]\n",
      " [345. 835. 147.]\n",
      " [495. 850.  98.]]\n",
      "矩阵A乘矩阵B：\n",
      "[[273]\n",
      " [339]\n",
      " [217]]\n",
      "矩阵A与矩阵c每个元素相乘：\n",
      "[[400  16   4  49]\n",
      " [100 324  25  36]\n",
      " [ 25  49 256   9]]\n",
      "A矩阵的逆矩阵是：\n",
      "[[ 3. -2.]\n",
      " [-1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "A = np.array([[150,250,50],[250,500,100],[300,700,120],[450,850,80]])\n",
    "B = np.array([[180,350,60],[300,550,120],[350,850,150],[500,850,100]])\n",
    "C = A + B\n",
    "print(f'矩阵A 加 矩阵B: \\n {C}')\n",
    "D = 0.1*A + 0.9*B\n",
    "print(f'矩阵A数乘 加 矩阵B数乘：\\n {D}')\n",
    "A = np.array([[20,4,2,7],[10,18,5,6],[5,7,16,3]])\n",
    "b = np.array([[8],[10],[5],[9]]) # 等价与 np.matrix('8;10;5;9')\n",
    "C = np.dot(A,b)\n",
    "print(f'矩阵A乘矩阵B：\\n{C}')\n",
    "A = np.array([[20,4,2,7],[10,18,5,6],[5,7,16,3]])\n",
    "c = np.array([[20,4,2,7],[10,18,5,6],[5,7,16,3]])\n",
    "D = A * c\n",
    "print(f'矩阵A与矩阵c每个元素相乘：\\n{D}')\n",
    "A = np.array([[1,2],[1,3]])\n",
    "X = np.linalg.inv(A)\n",
    "print(f'A矩阵的逆矩阵是：\\n{X}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 用最简行阶梯求逆矩阵\n",
    "1. ($A$,$I$)-->($I$,$A^{-1}$) 定理2.1推导出来的"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A的逆矩阵是：\n",
      " [[1. 1. 3.]\n",
      " [2. 1. 1.]\n",
      " [3. 2. 3.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,3,-2],[-3,-6,5],[1,1,-1]])\n",
    "X = np.linalg.inv(A)\n",
    "print(F\"A的逆矩阵是：\\n {X}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  3. -2.  1.  0.  0.]\n",
      " [-3. -6.  5.  0.  1.  0.]\n",
      " [ 1.  1. -1.  0.  0.  1.]]\n",
      "\n",
      " b的逆矩阵是：\n",
      " [[1. 0. 0. 1. 1. 3.]\n",
      " [0. 1. 0. 2. 1. 1.]\n",
      " [0. 0. 1. 3. 2. 3.]]\n",
      "排除左边的单位矩阵\n"
     ]
    }
   ],
   "source": [
    "b = np.eye(3)  # 生成3阶单位矩阵\n",
    "z = np.hstack((A,b))  # 水平合并A,b矩阵\n",
    "print(z)\n",
    "X = rref(z)  # 对b进行初等行变换\n",
    "print(F\"\\n b的逆矩阵是：\\n {X}\\n排除左边的单位矩阵\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.        , -1.        ,  0.        ,  1.        ,\n",
       "        -2.        ],\n",
       "       [ 0.        ,  1.        , -0.33333333,  0.        , -0.66666667,\n",
       "         1.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  1.        ,  1.        ,\n",
       "        -1.        ]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1,3,-2],[-3,-6,5],[-2,-3,3]])\n",
    "b = np.eye(3)\n",
    "z = np.hstack((A,b))\n",
    "rref(z)\n",
    "# 这个前面不是单位矩阵就没有逆矩阵了"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.25179981e+15, -2.25179981e+15,  2.25179981e+15],\n",
       "       [-7.50599938e+14, -7.50599938e+14,  7.50599938e+14],\n",
       "       [-2.25179981e+15, -2.25179981e+15,  2.25179981e+15]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1,3,-2],[-3,-6,5],[-2,-3,3]])\n",
    "np.linalg.inv(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LU分解\n",
    "1. $L^{-1}A$ = U  ==> $A = LU$\n",
    "2. 把增广矩阵分解成LU,下三角矩阵（Lower Triangle Matrix）、上三角矩阵（Upper Triangle Matrix）\n",
    "3. 前提方正、可逆\n",
    "> 参考：https://zhuanlan.zhihu.com/p/54943042\n",
    "\n",
    "问题：为什么LU分解速度比高斯消元法快？？？计算机没搞懂"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "增广矩阵C是：\n",
      " [[ 3  2 -2 -4]\n",
      " [ 3  3 -1 -5]\n",
      " [ 2  2 -1  4]]\n",
      "P是：\n",
      " [[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n",
      "L是：\n",
      " [[1.         0.         0.        ]\n",
      " [1.         1.         0.        ]\n",
      " [0.66666667 0.66666667 1.        ]]\n",
      " U是：\n",
      " [[ 3.          2.         -2.         -4.        ]\n",
      " [ 0.          1.          1.         -1.        ]\n",
      " [ 0.          0.         -0.33333333  7.33333333]]\n",
      "L*U是增广矩阵C：\n",
      " [[ 3.  2. -2. -4.]\n",
      " [ 3.  3. -1. -5.]\n",
      " [ 2.  2. -1.  4.]]\n"
     ]
    }
   ],
   "source": [
    "from scipy import linalg as la\n",
    "import numpy as np\n",
    "A = np.array([[3, 2,-2], [3, 3,-1],[2,2,-1]])\n",
    "b = np.array([[-4],[-5],[4]])\n",
    "C = np.hstack((A,b))\n",
    "print(f\"增广矩阵C是：\\n {C}\")  # C 就是 A\n",
    "P, L, U = la.lu(C)\n",
    "print(f\"P是：\\n {P}\")\n",
    "print(f\"L是：\\n {L}\\n U是：\\n {U}\")\n",
    "a = np.dot(L,U)\n",
    "print(f\"L*U是增广矩阵C：\\n {a}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 矩阵幂运算\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1,  2,  4,  8, 16], dtype=int32)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "A = np.array([2,2,2,2,2])\n",
    "B = np.array([0,1,2,3,4])\n",
    "np.power(A,B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 行列式运算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "直接算A的行列式值：\n",
      "-1660.0000000000011\n",
      "用LU=A分解算行列式值：\n",
      "1660.0000000000011\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg as la\n",
    "A = np.array([[-6,-7,6,7],[3,4,2,3],[-4,-2,0,6],[1,7,8,3]])\n",
    "# det(A) = det(L)det(U)、A = LU\n",
    "B = np.linalg.det(A)\n",
    "print(f\"直接算A的行列式值：\\n{B}\")\n",
    "P, L, U = la.lu(A)\n",
    "C = np.linalg.det(L)*np.linalg.det(U)\n",
    "print(f\"用LU=A分解算行列式值：\\n{C}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "273.067px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
